Thermodynamic anomalies in a lattice model of water 
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We investigate a lattice-fluid model of water, defined on a three-dimensional body centered cubic 
lattice. Model molecules possess a tetrahedral symmetry, with four equivalent bonding arms, aiming 
to mimic the formation of hydrogen bonds. The model is similar to the one proposed by Roberts 
and Debenedetti [J. Chem. Phys. 105, 658 (1996)], simplified in that no distinction between 
bond "donors" and "acceptors" is imposed. Bond formation depends both on orientation and local 
density. In the ground state, we show that two different ordered (ice) phases are allowed. At finite 
temperature, we analyze homogeneous phases only, working out phase diagram, response functions, 
the temperature of maximum density locus, and the Kauzmann line. We make use of a generalized 
first order approximation on a tetrahedral cluster. In the liquid phase, the model exhibits several 
anomalous properties observed in real water. In the low temperature region (supercooled liquid), 
there are evidences of a second critical point and, for some range of parameter values, this scenario 
is compatible with the existence of a reentrant spinodal. 

PACS numbers: 61.20.-p, 64.60.Cn, 64.60.My, 65.20.+W 



I. INTRODUCTION 

It is a well known fact that several thermodynamic 
properties of water exhibit some anomalous behaviori^. 
First of all, the heat capacity is unusually large and, at 
ordinary pressures, the solid phase (ice) is less dense than 
the liquid. Moreover, the liquid phase displays a temper- 
ature of maximum density at constant pressure, while 
both isothermal compressibility and isobaric heat capac- 
ity have a minimum as a function of temperature. Gener- 
ally speaking, the anomalous properties can be explained 
by the ability of water molecules to form hydrogen bonds, 
and by the peculiar features of such kind of bonds 4 ^. The 
same physics is thought to underly the unusual properties 
of water as a solvent for apolar compounds^, that is of 
the hydrophobic effect, whose importance in biophysics 
has been recognized in the latest year A Nevertheless, 
a comprehensive theory which explains all of these phe- 
nomena has not been developed yet. 

"Realistic" simulations of wate»2ii2iii*i£, based on 
more and more refined (but still phenomenological) in- 
teraction potentials, have reached quite a high level of 
accuracy in describing water thermodynamics. Neverthe- 
less, they are intrinsically limited by the large computa- 
tional effort required, which becomes still larger when it 
is necessary to determine multiple derivatives of the free 
energy, such as response functions. Moreover, due to the 
high level of microscopic detail, both of geometry and of 
interactions, often they do not make it easy to discrimi- 
nate what is essential to explain macroscopic properties. 
On the contrary, simplified models need simpler numer- 
ical calculations and, even if their quantitative accuracy 
is often poor, it is generally easier to trace connections 
between microscopic interactions and macroscopic prop- 
ertiea 13 ! 14 . 15 ! 16 ! 17 ! 18 . 19 . 20 ! 21 ! 22 ! 23 ! 24 . A simplified media- 
nism, proposed to account for the significant anomalies 
of water is the following one (see for instance Refs. 111 1251) . 



The formation of a hydrogen bond requires that two 
molecules assume certain relative orientations, staying at 
a distance larger than the one needed to minimize Van 
der Waals energy. This fact gives rise to a competition 
between the two kinds of interaction. Optimizing Van 
der Waals interaction allows higher density and higher 
orientational entropy, but yields a weaker binding energy, 
whereas, optimizing hydrogen bonding requires a lower 
density and a lower orientational entropy, but gives rise 
to a stronger binding energy. Therefore, at low enough 
temperature, local density and entropy fluctuations may 
become positively correlated, thus rationalizing a change 
of sign of the thermal expansion coefficient, that is a den- 
sity maximum. Such a simple mechanism has been im- 
plemented by different models, both on 22 i 23 i 26 i 27 i 28 anc j 
off-lattice^, in 3 22 i 23 as well as 2 dimensiona 24 : 26 ' 27 : 28 . 

One of them is the 3-dimensional model proposed by 
Roberts and Debenedetti (RD) 23-29 , defined on the body 
centered cubic lattice. Model molecules possess four 
bonding arms (two donors and two acceptors) arranged 
in a tetrahedral symmetry. Working on a lattice, one has 
to resort to a trick to describe hydrogen bond weaken- 
ing, when the two participating molecules are too close 
to each other. Such a trick is defined as follows. The 
energy of any formed bond is increased of some fraction 
(weakened bond) by the presence of a third molecule on 
a site close to the bond. Let us notice that the model has 
the same bonding properties as the early model proposed 
by Bell 1 ^, but the weakening criterion is different. The 
RD model is quite appealing in that it has been shown 
to predict some of real water thermodynamic anomalies, 
such as the temperature of maximum density, also show- 
ing evidence of a liquid-liquid phase separation in the su- 
percooled region, and of a second critical point. In view 
of investigations on mixtures of water with other chemi- 
cal species, as is the case, for instance, in most biological 
processes, it would be desirable to obtain an even simpler 
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FIG. 1: Two conventional cells of the body centered cubic 
lattice: A, B, C, D denote 4 interpenetrating face centered cu- 
bic sublattices. 




FIG. 2: Two model molecules forming a H bond. The lower 
molecule is in the i — 1 configuration, the upper one is in the 
i = 2 configuration. 



model, capable of capturing the same essential features. 
As it has been pointed out by other authors2431, the dis- 
tinction between donors and acceptors is likely to be not 
so crucial to describe the physics of hydrogen bonding. 
Therefore, a simplified version of the RD model, without 
a distinction between hydrogen bond donors and accep- 
tors, might be a good compromise between simplicity and 
accuracy. 

In this paper we investigate such a model, with a 
twofold purpose. As mentioned above, we are meant to 
explore the possibility of obtaining a simpler model with 
the same underlying physical mechanism, and with qual- 
itatively the same macroscopic properties. In addition, 
we are interested in providing a more detailed analysis of 
the effect of the weakening parameter, which turns out 
to be extremely relevant to determine the phase diagram, 
mainly in the supercooled liquid region. The paper is or- 
ganized as follows. In Sec. II we define the model and 
analyze its ground state. In Sec. Ill we introduce the 
first order approximation in a cluster variational formu- 
lation (cluster-site approximation), which we employ for 
the analysis. Sec. IV describes the results and Sec. V is 
devoted to some concluding remarks. An Appendix re- 
ports the calculation of density response functions and 
spinodals for the liquid phase. 



II. THE MODEL AND THE GROUND STATE 

Let us first introduce the model. Molecules are placed 
on the sites of a body centered cubic lattice, whose struc- 
ture is sketched in Fig. ^ A site may be empty or occu- 
pied by a water molecule. An attractive potential energy 
— e < is assigned to any pair of nearest neighbor (NN) 



occupied sites. This is the ordinary Van der Waals con- 
tribution. In the RD model, water molecules possess four 
arms that can form hydrogen (H) bonds (two donors and 
two acceptors), arranged in a tetrahedral symmetry, so 
that they can point towards 4 out of 8 NNs of a given 
site. We assume for simplicity that donors and accep- 
tors are undistinguishablc, that is a H bond is formed 
whenever two NN molecules have a bonding arm point- 
ing to each other, yielding an energy —77 < 0. Without 
such a distinction, it turns out that a water molecule 
has only 2 different configurations in which it can form 
H bonds (see Fig. |2J). We assume that w more configu- 
rations are allowed, in which the molecule cannot form 
bonds. The w parameter is related to the bond-breaking 
entropy. Moreover, to account for the fact that H bonds 
are most favorably formed when water molecules are lo- 
cated at a certain distance, larger than the optimal Van 
der Waals distance, the RD model assigns an energy in- 
crease r]c/6, with c € [0, 1], for each of the 6 sites closest 
to the bond occupied by a water molecule (i.e., 3 out of 6 
second neighbors of each participating molecule). A bond 
surrounded by all 6 water molecules is "fully weakened" 
and contributes an energy —77(1 — c). 

The hamiltonian of the system can be written as a sum 
over irregular tetrahedra, whose vertices lie on 4 different 
face-centered cubic sublattices, shown in Fig. ^ O ne of 
such tetrahedra is shown in Fig. E^a). We have 

^ = g y ] tiiaipi^isi (1) 
(a,/3,7,<5> 

where TLijki is a contribution which will be referred to as 
tetrahedron hamiltonian, and the subscripts i a ,ip,i 7 ,is 
label site configurations for the 4 vertices a, j3, 7, 6, re- 
spectively. Possible configurations are: Empty site (i = 
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FIG. 3: (a) Basic cluster (irregular tetrahedron): A,B,C,D 
denote sites in the 4 corresponding sublattices. AB, BC, CD, 
and DA are NN pairs; AC and BD are second neighbor pairs, 
(b) Husimi tree structure corresponding to the generalized 
first order approximation on the tetrahedron. 



0), site occupied by a molecule in one of 2 bonding orien- 
tations (i — 1,2; see Fig.|2l or in one of w non-bonding 
configurations (i — 3). Assuming that (j, k), (k,l), 

and (/, i) refer to NN pair configurations, the tetrahedron 
hamiltonian reads 
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ijkl 
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where rij is an occupation variable, defined as n, = 
for i = (empty site) and m = 1 otherwise (occu- 
pied site), while hij = 1 if the pair configuration 
forms a H bond, and hij = otherwise. Let us also as- 
sume that k, I (in this order) denote configurations of 
sites placed on, say, A, B, C, D sublattices respectively. If 
A,B,C,D are defined as in Fig.^ we can define hij = 1 
if i = 1 and j = 2, and hij = otherwise. With the 
above assumptions, the tetrahedron hamiltonian is inde- 
pendent of the orientation, that is of the arrangement of 
A, B, C, D on its vertices. This is a nice effect of removing 
the donor-acceptor distinction, which considerably sim- 
plifies the whole analysis. Let us notice that both Van 
der Waals (— erijnj) and H bond energies (— which 
are 2-body terms, are split among 6 tetrahedra, whence 
the 1/6 prefactor in Eq. (|TJ . On the contrary the 3-body 



weakening terms (ijchijUk/d) are split between 2 tetra- 
hedra, thus the 1/6 factor is absorbed in the prefactor, 
while a 1/2 factor is left in the tetrahedron hamiltonian. 
Let us denote the tetrahedron configuration probability 
by pijki , with the same convention about the subscript or- 
der, and assume that the probability distribution is equal 
for every tetrahedron. Taking into account that there are 
6 tetrahedra per site, we can write the following expres- 
sion for the internal energy per site of an infinite lattice 



3 3 3 3 

i=0 1=0 fc=0 1=0 



(3) 



The multiplicity for the tetrahedron configuration 



(i,j,k,l) is given by WiWjWkWi, where Wi = 



for 



i = 3 (non-bonding configuration) and Wi = 1 otherwise 
(bonding configuration or vacancy). 

Let us now have a look at the ground state of the 
model. In order to do so, let us investigate the zero 
temperature grand-canonical free energy uj° = u — /ip 
(/i being the chemical potential and p the density, i.e., 
the average site occupation probability), which can be 
formally written in the same way as the internal energy u 
of Eq. J2J), replacing the tetrahedron hamiltonian 'Hijki 
by 



H 



ijkl 



H 



ijkl fl>~ 



rii + rij + n k + ni 



(4) 



First of all, we have an infinitely dilute "gas" phase (G) 
with zero density and zero free energy. Secondly, we can 
conceive an ordered "open" ice phase (0) with two fully 
occupied sublattices, say A and B, while the other sub- 
lattices, C and D, are empty (p — 1/2). In this con- 
figuration all molecules are fully bonded, i.e., there are 
2 bonds per molecule, and no bond is weakened, so that 
the free energy turns out to be 



LOn = -e 



V - M/2. 



(5) 



Another possibility is the "closed" ice phase (C), in which 
all sites are occupied (p = 1), the maximum number of 
H bonds is formed (for instance, all AB and CD pairs are 
bonded), but all bonds are fully weakened. The resulting 
free energy is 

uj° c = -4e - 277(1 -c)-fi. (6) 

It is easy to show that the G phase is thermodynamically 
favored (0 < ujq and < w£) for fx < pqo, where 

MGO - -2c - 277, (7) 

the O phase is favored (wq < and ujq < ujq) for pqq < 
p, < poc, where 

-6e - 277(1 - 2c), 



Moc 



< UJ°r 



(8) 
for 



and the C phase is favored (ujq < and ui, 
M > /^OC- The O phase has actually a stability region 
i.e., pqo < floe, provided 



c > e/77. 



(9) 
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We shall always work in the latter regime, which allows 
to reproduce two different forms of ice. Even if in the 
following we shall not deal with ordered phases at finite 
temperature, the latter choice should be the most reason- 
able one, in order to describe real water properties. We 
have considered also the possibility of different structured 
phases, respectively with 1 or 3 occupied sublattices, but 
they never turn out to be stable, in the physical range of 
parameter values. 

III. FIRST ORDER APPROXIMATION 

We perform the finite temperature analysis by means 
of a generalized first order approximation on a tetrahe- 
dron cluster. Let us introduce the approximation in the 
framework of the cluster variation method, an improved 
mean-field theory which in principle can take into ac- 
count correlations at arbitrarily large, though finite, dis- 
tances. In Kikuchi's original work^I, an approximate en- 
tropy expression was obtained by heuristic arguments, 
while, in more recent and rigorous formulations^, the 
approximation is shown to be equivalent to a trunca- 
tion of a cluster cumulant expansion of the entropy. The 
approximation is expected to work, because of a rapid 
decreasing of the cumulant magnitude, upon increasing 
the cluster size, namely when the latter becomes larger 
than the correlation length of the system^ 3 -. A particular 
approximation is defined by the largest clusters left in 
the truncated expansion, usually denoted as basic clus- 
ters. One obtains a free energy functional in the cluster 
probability distributions, to be minimized, according to 
the variational principle of statistical mechanics. 

For our model we choose a number of irregular tetra- 
hedra as basic clusters, namely 4 out of 24 tetrahcdra 
sharing a given site, as sketched in Fig. [21 This choice 
actually turns out to coincide with the (generalized) first 
order approximation (on the tetrahedron cluster), which 
is also equivalent to an exact calculation on a Husimi 
latticed, whose (tetrahedral) building blocks are just ar- 
ranged as in Fig. ISIb). Such an approximation has not 
only the advantage of high simplicity, due to the fact 
that the only clusters retained in the expansion are ba- 
sic clusters and single sites (it is sometimes referred to 
as cluster-site approximation 35 ), but also of a relative 
accuracy, which has been recognized for different mod- 
els, even with orientation dependent interactions^. Let 
us notice that the internal energy is treated exactly, be- 
cause the range of interactions does not exceed the basic 
cluster size. The grand canonical free energy per site 
uj = u — Ts — up (T being the temperature and s the 
entropy per site) can be written as 

u 3 3 3 3 \ 

- = w * Y w i Y Wk Y w iPW I —^r- + ^Pm I 
1=0 j=Q k=0 i=o \ / 

3 

-3^ Wjpi \npi , (10) 

i=0 



where pi is the probability of the site configuration i 
(temperature is expressed in energy units, whence en- 
tropy in natural units). In this paper we focus on liq- 
uid, i.e., homogeneous state, hence we assume that all 
sites have the same configuration probability distribu- 
tion. The latter can then be obtained as a marginal of 
the tetrahedron distribution pijki, by the following sym- 
metrized expression 

3 3 3 

E\ -> \ - Pijkl + Plijk + Pklij + Pjkli 
W J 2^ Wk l^ Wl 4 ■ 

j=0 k=0 1=0 

(11) 

The free energy turns out to be a function of the only 
tetrahedron probability distribution, taken as variational 
parameter. The minimization with respect to such pa- 
rameter, with the normalization constraint 

3 3 3 3 

Y m Y w i Y Wk Y w iPvM = i, (12) 

i=0 j=0 k=0 1=0 

can be performed by the Lagrange multiplier method, 
yielding the equations 

Pijkl = C X e-^ kl,T {p l p J PkPif i , (13) 

where £, related to the Lagrange multiplier, can be com- 
puted by imposing the constraint Eq. (|12fl as 

3 3 3 3 

6 = Y m Y w i Y Wk Y w i e ~ n * jki,T {piPjPkPif /A ■ 

i=0 j=0 k=0 1=0 

Eq. (|13fl is in a fixed point form, and can be solved numer- 
ically by simple iteration (natural iteration method^). 
For the cluster-site approximation, the numerical proce- 
dure can be proved to reduce the free energy at each iter- 
ation 3 ^ 3 ^, and therefore to converge to local minima. Let 
us notice that the symmetrized marginalization Eq. i|ll|) 
imposes implicitly a homogeneity constraint, useful to 
investigate the metastable liquid properties. In fact it is 
easy to see that Eq. (|13fl , due to invariance of TLijki under 
cycle permutation of the subscripts [see Eqs. 10 and (@J], 
determines a tetrahedron distribution Piju with the same 
property, giving rise to four equal site marginals. There- 
fore, we must be aware that the stationary point we find 
may be unstable with respect to a translational symme- 
try breaking. Anyway, from the solution of Eq. (|13jl one 
obtains a tetrahedron probability distribution {pijkl}, 
whence one can compute the thermal average of every 
observable, the internal energy by Eq. @ and the free 
energy by Eq. (fTO)) . The latter can be also related to the 
normalization constant as 

w = -Tln£, (15) 

whence £ can be viewed as an effective (single site) grand 
canonical partition function. We shall make use of this 
property in the following. 
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IV. THERMODYNAMIC PROPERTIES 

In order to investigate the model properties, let us fix 
a set of parameters. First of all, let us take e/r] = 0.25. 
This value is equal to the one employed for a previous 
mean field analysis^, and similar to the one chosen by 
RD for the original model 2 ^. This choice accounts for the 
greater binding energy of hydrogen bonds with respect 
to Van der Waals interactions, and will be kept fixed 
throughout the present analysis. As far as the multiplic- 
ity of non-bonding water configurations is concerned, we 
set w — 20, to mimic the high directionality of hydrogen 
bonds. From the phase diagram analysis, it turns out 
that it is necessary to set this parameter large enough 
to let anomalous properties appear, but further increase 
does not change qualitatively the phase behavior. There- 
fore, also the latter parameter will be held fixed in the 
following. On the contrary, we shall investigate in detail 
the effect of changing the weakening parameter c, which 
is actually the crucial one for the present model, in the 
regime c > 0.25, according to Eq. 10. 

A. Phase diagrams 

Let us report in Fig. 0] a sequence of temperature- 
pressure (left column) and density-temperature (right 
column) phase diagrams, for different c values. Imposing 
homogeneity, our analysis includes both thermodynami- 
cally stable and metastable (supercooled) phases, even if 
stability is not investigated. Let us notice that pressure 
can be determined as P = —u, due to the fact that the 
free energy has been defined as a grand-canonical poten- 
tial. We have assumed volume per site equal to 1, i.e., 
pressure is expressed in energy units. 

For c = 0.3 we have essentially an ordinary liquid- 
vapor coexistence line, terminating at a critical point, 
even if an anomalous temperature dependence of the liq- 
uid density (without a maximum) can be observed. In 
the very low temperature region, an intermediate density 
liquid phase appear, giving rise to a triple point. This 
phase is actually an unphysical solution of our equations, 
in that the entropy turns out to be negative. In this re- 
gion a crystalline phase is likely to be stable. 

For c = 0.4 the phase diagram undergoes a dramatic 
change. The intermediate liquid phase region becomes 
a unique region with the ordinary liquid, and the triple 
point is replaced by a second critical point, where the low- 
high density liquid coexistence terminates. The second 
critical point still lies in a negative entropy region. The 
density of the liquid coexisting with the vapor displays a 
maximum as a function of temperature. 

For c = 0.5 the density maximum becomes more pro- 
nounced, but there is no topological change in the phase 
diagram. 

For c = 0.6 a new coexistence line appears between 
the high density and the low density liquid. Such a tran- 
sition line (probably metastable) is negatively sloped in 
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FIG. 4: Pressure (P/rj) vs temperature (T/rj) phase diagrams 
(left column) and temperature vs density (p) phase diagrams 
(right column) for e/r/ = 0.25 and w = 20. From top to 
bottom c = 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, respectively. Thick solid 
lines denote (first order) phase transitions (left) and delimit 
coexistence regions (right). A thin dashed line (right) corre- 
sponds to three-phase coexistence (triple point). 



the temperature-pressure phase diagram, and terminates 
at a third critical point. In this scenario, the ordinary liq- 
uid phase stability is delimited by a reentrant spinodal, 
as will be verified in the following. A similar scenario has 
been already observed in a two-dimensional models. 

For c = 0.7 the two "metastable" critical points gets 
closer, and finally, for c = 0.8, they merge, giving rise 
to a unique high-low density liquid coexistence, which 
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terminates at a triple point. The topology of the phase 
diagram is again similar to the one obtained for low c 
values, but in this case there exists a region where the low 
density liquid has positive entropy. From a microscopical 
point of view, such a liquid phase is characterized by a 
high probability of bonding configurations, i.e., it is a 
highly hydrogen bonded phase. This feature corresponds 
to a lower density, as previously pointed out. 

The role of the weakening parameter is well charac- 
terized by the sequence of temperature-pressure phase 
diagrams. Actually, the general trend is that, upon in- 
creasing c, the stability region of the low density liquid 
(having a low number of weakening molecules) gets larger 
and larger. In the next part of this work we focus on two 
particular parameter choices (c = 0.4, 0.6), as representa- 
tives of a range of c values in which the model results are 
qualitatively consistent with the anomalies observed for 
real liquid water in ordinary temperature and pressure 
conditions. The general analysis reported above is com- 
pleted by studying the locus of density extrema (max- 
ima and minima) as a function of temperature, the liq- 
uid phase spinodals, and the Kauzmann line, where the 
liquid entropy vanishes (ideal glass transition). 



B. TED locus, spinodal, and Kauzmann line 

One of the thermodynamic anomalies of the present 
model is the temperature of maximum density (TMD) 
along isobars for the liquid phase. Nevertheless, for some 
pressure range, it is possible to observe also a temper- 
ature of minimum density, therefore we shall generally 
speak about a temperature of extremum density (TED). 
Joining temperatures of maximum (or extremum) den- 
sity at different pressures defines the so called TMD (or 
TED) locus. At ordinary pressure, the TMD locus is a 
negatively sloped line in the T-P phase diagram of real 
water. In principle, we could determine the TED lo- 
cus numerically, by adjusting the chemical potential in 
order to fix the pressure and then imposing that the (iso- 
baric) thermal expansion coefficient vanishes. Actually 
we have performed a different calculation, based on the 
effective partition function l|15p. rewritten as a function 
of only two variational parameters, namely, the density 
p and the fraction <fi of bonding molecules. Details about 
this calculation, which allows to determine density re- 
sponse functions and spinodals as well, are given in the 
Appendix. The limit of stability of the liquid phase (spin- 
odal) is the locus in which the metastable liquid ceases 
to be a minimum of the free energy, and becomes a sad- 
dle point. Actually, in the homogeneity hypothesis, we 
add a constraint to the free energy, thus neglecting sta- 
bility loss with respect to symmetry broken (ice) phases. 
Aware of this, spinodals can be obtained by imposing 
that the hessian determinant of the homogeneous free 
energy (or partition function) vanishes, as shown in the 
Appendix. Let us notice that, in this case, it is not pos- 
sible to work out the result making use of the natural 
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FIG. 5: Pressure (P/rj) vs temperature (T/rj) phase diagram 
(top panel) and temperature vs density (p) phase diagram 
(bottom panel) for e/77 = 0.25, w = 20, c = 0.4. Thick solid 
lines denote (first order) phase transitions (top) and delimit 
coexistence regions (bottom). Thin (solid, dotted, and dash- 
dotted) lines denote spinodals, TED locus, and Kauzmann 
line, respectively. 



iteration Eqs. because the attraction basin of the 

liquid phase gets smaller and smaller, and vanishes at 
the spinodal. On the contrary, the locus at which the 
liquid phase entropy vanishes (Kauzmann line), can be 
easily determined numerically. 

The results are shown in Figs. and for c = 0.4 
and c = 0.6, respectively. In both cases, we find a den- 
sity maximum as a function of temperature for liquid 
coexisting with vapor (and at constant pressure as well), 
and the TMD slightly decreases upon increasing pressure. 
Another common feature is the presence of a high-low 
density liquid coexistence, and of a corresponding critical 
point. For c = 0.4 the critical point lies in the negative 
entropy region, beneath the Kauzmann line, while for 
c = 0.6 it lies in the positive entropy region. The Kauz- 
mann line displays a cusp (towards low temperature) in 
the vicinity of the critical point. The low approximation 
level of our analysis does not allow to take all of this 
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FIG. 6: Pressure (P/rj) vs temperature (T/rj) phase diagram 
(top panel) and temperature vs density (p) phase diagram 
(bottom panel) for e/rj — 0.25, w — 20, c = 0.6. Thick solid 
lines denote (first order) phase transitions (top) and delimit 
coexistence regions (bottom). A thin dashed line (bottom) 
corresponds to three-phase coexistence (triple point). Thin 
(solid, dotted, and dash-dotted) lines denote spinodals, TED 
locus, and Kauzmann line, respectively. 
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FIG. 7: Pressure (P/rj) vs temperature (T/rj) phase diagram 
for e/rj = 0.25, w = 20, c = 0.5. A thick solid lines denotes 
the liquid-vapor transition. The thin (solid and dotted) lines 
denote the liquid spinodal and the TMD locus, respectively. 



c = 0.6 the maximum occurs at a far higher pressure than 
for c = 0.4. After the maximum, upon decreasing tem- 
perature, the TMD locus becomes a temperature of min- 
imum density locus, ending in the second critical point. 
After that, the line continues (again as a TMD locus) in 
the metastability region of the high density liquid with 
respect to the low density liquid. Upon decreasing pres- 
sure, we find the most relevant differences between the 
two cases. For c = 0.4, the TMD reaches a maximum and 
then decreases again, as we move in the negative pressure 
region. For c = 0.6, the TMD always increases. In both 
cases the TMD locus terminates against the liquid-vapor 
spinodal. In the T-P diagram, the meeting point corre- 
sponds to a minimum of the spinodal line, as required by 
the thermodynamic consistency arguments of Speedy and 
Debenedetti^2iiMiii2i^. Such a minimum is placed at a 
much higher temperature for c = 0.6 than for c = 0.4. 
As far as the liquid- vapor spinodal is concerned, another 
important difference is observed, namely, for c = 0.6 the 
spinodal reenters the positive pressure region, ending in 
the above mentioned "third" critical point. Such a reen- 
trance gives rise to a divergence in the response functions, 
measured at constant pressure, as we shall see below. 
This is not the case for c = 0.4. As one could expect, it 
is possible to verify that the spinodal reenters the posi- 
tive pressure region, if and only if the third critical point 
exists. The boundary between the two regimes is found 
to be around (actually a little greater than) c = 0.5. A 
region of the T-P phase diagram for c = 0.5, showing the 
reentrance of the TMD locus, is reported in Fig. [7| 



information as reliable, but let us notice that a similar 
scenario has been predicted also by a statistical analy- 
sis of the potential energy landscape of simulated wa- 
ter, performed by Sciortino and coworkers on the basis 
of the inherent structure theory and of some simplifying 
assumptions 38 . As far as the TMD locus is concerned, 
we observe a peculiar maximum in the T-P plane. For 



C. Response functions 

Let us now investigate the density response functions 
and the specific heat of the liquid at constant pressure 
P/rj = 0.05, roughly corresponding to 1/10 of the liquid- 
vapor critical pressure. Also for this calculations, dc- 
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0.28 0.30 0.32 0.34 0.36 0.38 
temperature 



FIG. 8: Response functions at constant pressure (P/rj = 
0.05) for the liquid phase as a function of temperature (T/rf), 
for e/r/ = 0.25 and w = 20. From top to bottom we show the 
isobaric thermal expansion coefficient (rjotp), the isothermal 
compressibility (t)Kt), and the specific heat (cp). Numerals 
beside each plot denote c values. 



tails are left to the Appendix. We find anomalous be- 
havior, similar to that of real liquid water. The first 
response function we consider is the thermal expansion 
coefficient up = (—d In p/dT)p, which, from statistical 
mechanics, is known to be proportional to the entropy- 
volume cross-correlation. For ordinary fluids, ap is al- 
ways positive, i.e., the local entropy and the local spe- 
cific volume are positively correlated. On the contrary, 
for our model ap (Fig. |SJ top panel) is anomalous. As 
temperature is lowered, the expansion coefficient van- 
ishes (at the TMD), and then becomes negative. For 
c = 0.4, we have observed just a shallow minimum (not 
shown) around T/rj — 0.17. The minimum is more pro- 
nounced for c = 0.5, which, as mentioned above, is still 
in a regime where the third critical point does not ex- 
ists. Finally, for c = 0.6, the minimum becomes a di- 



vergence, due to the spinodal reentrance in the positive 
pressure region. Of course, the divergent behavior can 
be observed only for pressure values less than the third 
critical point pressure. The trend of the isothermal com- 
pressibility kt — (din p/8P)t is also anomalous (Fig. [SJ 
middle panel). For a typical liquid, kt decreases as one 
lowers temperature, because it is proportional to density 
fluctuations, whose magnitude decreases, upon decreas- 
ing temperature. On the contrary, we can observe that 
kt, once reached a minimum, begins to increase upon 
decreasing temperature. Only a broad maximum is ob- 
served for c = 0.4; the maximum becomes sharper for 
c = 0.5, and finally becomes a divergence for c = 0.6. 
The constant pressure specific heat cp = (~T d 2 p / dT 2 ) p 
(Fig. [SJ bottom panel) displays a completely analogous 
behavior, with the minimum occurring at a higher tem- 
perature. Qualitatively similar thermodynamic anoma- 
lies arc observed in real liquid water, even if the possibil- 
ity of divergent-like behavior in the supercooled regime 
is no longer believed to be realistic 3 . 



D. A comparison with the RD model 

Let us finally report the results of a comparison with 
the RD model, of which, as previously mentioned, the 
present model is a simplified version. In particular, let 
us consider the results of a Monte Carlo simulation 29 , 
in which parameters were chosen to push the low-high 
density liquid critical point to a temperature of the same 
order of magnitude of the ordinary liquid-vapor critical 
point. In this work, evidence was given that the model 
actually predicts liquid-liquid coexistence, and that the 
latter is not an artifact of approximations. The model 
parameters are e/rj = 0.2, c = 0.8, while the total num- 
ber of molecule orientations is q = 108. In order to a 
comparison, the latter parameter is to be renormalizcd, 
in that RD model molecules, possessing bonding arms of 
two different kinds (donors and acceptors), allow 12 dif- 
ferent bonding configurations, while the present model 
allows only 2. We have performed the renormalization as 
follows 



obtaining w — 16. The results are reported in Fig.|5| The 
topology of the (temperature-density) phase diagram is 
qualitatively similar, and also the quantitative agreement 
of critical temperatures is good. Only the density of the 
liquid coexisting with the vapor displays a significant dis- 
crepancy. The reentrance in the liquid-liquid coexistence 
region at low temperature is reproduced, even if RD con- 
jecture that a lower critical point exists, mainly on the 
basis of first order approximation results^, whereas we 
find coexistence down to zero temperature. Anyway, as 
shown before, the presence or absence of a lower critical 
point may be induced by small c value variations. 
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FIG. 9: Temperature (T/rj) vs. density (p) phase diagram 
for e/?7 = 0.2, w — 16, c = 0.8 (solid lines), compared to 
Monte Carlo results from Ref . [2t] (scatters) . Dotted lines are 
eyeguides. 



V. DISCUSSION AND CONCLUSIONS 

In this paper we have investigated a lattice fluid model 
with water-like features, by a generalized first order ap- 
proximation. As mentioned in the Introduction, the 
mechanism by which our model describes water anoma- 
lies is essentially based on the competition between the 
Van der Waals isotropic interaction and the highly di- 
rectional H bonding interaction, and on the difference 
between the respective optimal interaction distances. In 
the lattice framework, the latter is taken into account 
by means of the weakening trick, proposed by RD 23 . 
We have actually simplified the original RD model by 
neglecting the distinction between H bond donors and 
acceptors. In spite of this, the model turns to repro- 
duce qualitatively several thermodynamic anomalies of 
real water at constant pressure: a maximum of density, 
and a minimum of isothermal compressibility and specific 
heat. This fact confirms the idea that the distinction be- 
tween donors and acceptors is not crucial to the physics 
of H bonding 2 ^. By a suitable parameter scaling, we 
have also compared some results of the simplified model 
with those of the original model obtained by Monte Carlo 
simulations, and we have verified that the simplification 
and the approximation seem to preserve the qualitative 
features of the phase diagram, though in a case of quite 
simple topology. On the simplified model, it is possi- 
ble to perform a generalized first order approximation on 
a tetrahedron cluster, without introducing further "ad 
hoc" approximations, as it was necessary in the original 
paper by RD 23 . Moreover, with respect to Monte Carlo 
simulations, the first order approximation has not only 
the advantage of a much lower computational effort, but 
also allows a well defined extrapolation of the equation of 
state in the supercooled regime, without the need of cri- 
teria to prevent the simulation dynamics from falling into 
crystalline states. In fact, analysis of the model predic- 



tions in the metastable phase region is extremely inter- 
esting to rationalize observed thermodynamic anomalies. 

As previously shown, there is evidence of a second 
(metastable) critical point, which terminates a line of 
first order transitions between two liquid phases at dif- 
ferent densities. For c < 0.5, the second critical point lies 
at a temperature lower than the Kauzmann temperature, 
at which the configurational entropy vanishes. Moreover, 
the Kauzmann line displays a reentrance in the vicinity 
of the second critical point. All of these features turn out 
to be in a remarkably good agreement with the simpli- 
fied statistical analysis of the potential energy landscape 
of simulated water, recently performed by Sciortino and 
coworkers 3 ^. 

We have also computed the temperature of extremum 
density locus and the liquid spinodals, in the framework 
of our approximation. In the ordinary temperature and 
pressure region, the TED locus is a negatively sloped line 
(in the T-P diagram) corresponding to a density maxi- 
mum, as observed in experiments. The locus displays a 
pressure maximum, and, for higher pressure values, the 
liquid becomes normal. After the maximum, upon de- 
creasing temperature, the TED line denotes a minimum 
density, which is a peculiar feature of our model. The 
TED line crosses the Kauzmann line very close to its 
reentrance, as recently predicted by Speedj^. On the 
other side, the TED line is reentrant for low c values and 
it is not for higher c values. In both cases it terminates at 
a minimum of the spinodal in the T-P plane, as required 
by thermodynamic consistency^. Nevertheless, in the 
former case the spinodal does not reenter the positive 
pressure half-plane, while in the latter case it does. In 
this case, response functions exhibit a divergent behavior. 
Such results are interesting in that the reentrant spinodal 
was one of the conjectures invoked to explain thermody- 
namic anomalies of wates 3 ^^. The fact, that the same 
model with different parameter values may predict or not 
a reentrant spinodal, suggests that in real water there is 
probably a subtle balance of interactions, which does not 
allow to discriminate easily between the two possible sce- 
narios. This fact has been previously pointed out by phe- 
nomenological models^, but, to our knowledge, not yet 
by microscopical models. The most recent and accurate 
molecular dynamics simulations suggest a scenario with a 
non-reentrant spinodal and a reentrant ("nose shaped") 
TMD locus 3 -. The present model accounts for such a sce- 
nario for c « 0.5. In this case, a slight reentrance of the 
spinodal, induced by the vicinity of the TMD locus, can 
be observed. A similar "intermediate" situation is con- 
sistent with the results of an analytical equation of state 
derived by Truskett et al A We have remarked this fea- 
ture in that it may justify the early experimental results 
which supported the reentrant spinodal conjecture^. 

As a conclusion, let us remind that the model in the 
present treatment is not able to provide microscopic 
structural details as simulations do, but its most appeal- 
ing feature is simplicity. In spite of how little we put in 
the model, the latter turns out to give a qualitatively cor- 
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rect description of the peculiar thermodynamics of wa- 
ter, and is consistent with predictions based on much 
more sophisticated models and simulations. Therefore, 
we suggest that it may be suitable as a starting point to 
investigate complex phenomena, such as the solvation of 
apolar solutes, in which the physics of hydrogen bonding 
plays a significant role. 



APPENDIX: SPINODALS AND RESPONSE 
FUNCTIONS 

In this appendix we report in detail the calculation of 
the spinodal lines and of the density response functions, 
i.e., the (isobaric) thermal expansion coefficient ap and 
the isothermal compressibility kt- The (isobaric) heat 
capacity cp is determined by a numerical derivative of the 
chemical potential /i, which comes out in a natural way 
from the same calculation, as a function of temperature 
and pressure. 

Eq. I|15|l suggests that the normalization constant £ 
[Eq. 114(1 ] can be used as a variational form of the effec- 
tive (single site) grand canonical partition function, to be 
maximized with respect to the site probability distribu- 
tion pi . We can also observe that the homogeneous liquid 
phase is isotropic, whence the 2 bonding configurations 
on each site must have the same probability 



Pi =P2- 



(A.l) 



With the above assumption, we can replace the fourfold 
summation over the configurations i,j,k,l of a tetrahe- 
dron with a double summation over the number n = 
0, . . . , 4 of occupied sites and the number k = 0, . . . , n of 
molecules in a bonding configuration. There are ( ) (?) 
ways of arranging molecules with the above require- 
ments, and the multiplicity of the resulting configuration 
(k bonding molecules and n — k non bonding ones) is 
2 k w n ~ k , whence we obtain 



^ = EE 

n=0 fc=0 



k n „n—h / 4- nkri — k 



2 K w 



[Po P1P3 



3/4 



where x n ^k are linear combinations of the Boltzmann 

weights e~ nijk '/ T , appearing in Eq. ifH3|) . They are ob- 
tained as described in Tab. [I] In order to avoid imposing 
the normalization constraint, it is convenient to rewrite 
the site probabilities as a function of the density p (oc- 
cupation probability) and the fraction (f> of molecules in 
a bonding configuration. It is easy to obtain 



Pa 
Pi 

P3 



1 



P 



(A.3) 
(A.4) 

(A.5) 







2 
2 1 
2 2 



3 

3 1 

3 2 

3 3 



4 

4 1 

4 2 

4 3 

4 4 



re/4 



3 ^ 3 C 

§ + ¥ IT 

3- + §^ /T 



{| + I[I + ^ (W2)/T ]} 



1 _|_ l e v(l-c/2)/T 



S 4,/T 
2 4,/T 



1 -I- | e l(l-<0/T 

1 + e e v(i-c)/T + x e a^(i- c )/T| 



TABLE I: x n ,k coefficients. The way they have been com- 
puted is explained hereafter. The chemical potential is taken 
into account by the z n ^ A prefactor, where z = e M//T is the 
fugacity. For n = 0, 1 no interaction is possible, whence only 
the fugacity term is present. For n — 2 there are 6 possible 
arrangements of molecules on the sites; in 2 of them (1/3) the 
2 molecules are not NNs and there is no interaction; in the 
remaining 4 cases (2/3) the 2 molecules are NNs and interact 
with the Van der Waals energy (whence the factor e e ^ T ) . Only 
when k = 2 a H bond can occur (whence the factor e v ^ T ), in 
1 out of 4 possible configurations of 2 bonding molecules, i.e., 
when the molecules point an arm to each other. The H bond 
is not weakened because the neighbor sites are empty. For 
n = 3 all possible arrangements have 2 pairs of NN molecules, 
whence the factor e 2e ^ T . If k = 0, 1 no other interaction ex- 
ists. If k = 2 we can place the 2 bonding molecules on NN 
sites in 2/3 of the cases, and in 1/4 of their allowed configura- 
tion they form a H bond, which is half weakened, because one 
of the NN sites is occupied. With k = 3 bonding molecules, 
there are 2 3 = 8 possible configurations, 4 of which have 1 
formed bond. For n — 4 all sites are occupied and there are 
4 NN pairs, whence the factor e 4C//T . With k = 0, 1 bonding 
molecules, no other interaction is possible. With k = 2, there 
are 4/6 = 2/3 arrangements in which the bonding molecules 
are placed on NN sites and form a H bond in 1/4 of their 
configurations. The bond is fully weakened, because both the 
neighbor sites are occupied. With k = 3 bonding molecules, 
the explanation is equivalent to the n = 3 case. Finally, with 
k = 4 bonding molecules, it is necessary to enumerate all the 
2 4 = 16 configurations, 2 of which have no bond, 6 have 1 



bond, and 2 have 2 bonds. 



whence 



4 

E 

71 = 



E 

k=0 



J/ 4 



(A.6) 
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where 



and 



f n ,k(x) = x k (l-x) r 



(A.7) 



In order to impose thermodynamic equilibrium, we have 
to minimize £ with respect to the variational parameters 
p and <p. Let us notice that £ also depends on the temper- 
ature T and the chemical potential p (or equivalently the 
fugacity z — e^ T ), as appropriate in the grand-canonical 
ensemble. Such a dependence is hidden in the x n ^ coef- 
ficients (see Tab.HJ. We impose the necessary minimum 
condition, by setting the derivatives of £ with respect to 
p and <j) to zero. Moreover, we are interested in work- 
ing at fixed pressure P, therefore we have to impose a 
third equation u> = —P, where, from Eq. (|15fl . we have 
w = — Tlog£. We have thus to solve a system of three 
equations 



d(lnp) 

0OM) 

~e p ' T 



= 
= 

= o, 



(A.8) 

(A.9) 
(A.10) 



which is actually an implicit definition of the three func- 
tions p(T,P), <j>(T 7 P), and z(T,P). Having determined 
numerically the fugacity z(T,P), we immediately obtain 
the chemical potential p — Tlogz, whence the specific 
heat cp = (—Td 2 p/dT 2 )p by a numerical derivative. 

Let us now consider the three simultaneous Eqs. (|A.8(1 . 
(|A.9|I . and (jA.lOfl . As previously mentioned, the left 
hand sides are functions of p, <f>, z, T, P. If we replace 
p, (j), z by the previously determined functions of T, P, 
we obtain three functions that are identically equal to 0, 
for each value of T and P. Therefore, taking the par- 
tial derivatives with respect to X — T, P, we obtain the 
following linear system: 



A- 



/d(lnp)\ 




dX 




d(bx(t>) 




OX 




d(lnz) 




V dX J 





t 



where 
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9(lnp)9(ln (f>) 9(lnp)9(ln z) 
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d(ln<j)) 2 
9(ln0) 



d(]n<f>)d(]nz) 
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(A.12) 



d(P/T) 
dT 

d(P/T) 
dP 



P 



1 

T' 



(A.13) 
(A.14) 



The derivatives of £ can be easily determined from 
Eq. I|A.6|) . making use of 



d fn,k( x ) _ a{k-nx) 
d{hix) 



1-x 



d2 fn,k( x ) a 2 (k~nx) 2 -a{n-k)x 



d{hix) 



(1 



(A.15) 
(A.16) 



obtained from Eq. I|A.7(1 . a being a generic exponent, and 
of 



d{\nz) 



n 

~7 •En,ki 



(A.17) 



while dx n ^jdT can be obtained from Table ||| In or- 
der to determine the density response function, we solve 
Eq. (|A.11|) with respect to 9(ln p)/dX, making use of the 
Kramer's rule, yielding 



d(lnp) detB 



dX 



detA' 



(A.18) 



where B is obtained by replacing the right hand side of 
Eq. (jA.llfl in the first column of A. In this way the 
isobaric thermal expansion coefficient ap — —d(lnp)/dT 
and the isothermal compressibility kt — 9(ln p) / dP are 
expressed as a function of p,4>, z,T, P. As previously 
mentioned, p, 4>, z can be determined numerically as a 
function of T, P, for instance by solving the simultaneous 



Eqs. (JAjyi, (AJ5J, and ljAT0|) . As far as the TED locus 
is concerned, we have to add a fourth equation, imposing 
ap = 0. Finally, spinodal lines can be obtained as the 
locus in which the above response functions diverge, that 
is, detA = 0. Let us observe that, in the latter case, only 
three simultaneous equations are needed. In fact, if the 
grand canonical equilibrium conditions i|A.8|) and 1A.9I) 
are verified, the first two elements of the third line of A 
(which we may denote by ^31,^.32) vanish, therefore we 
can write det A = ^33 det A33, where A33 is the submatrix 
of A obtained by removing the third row and the third 
column. As a consequence, the equation det A = does 
not contain P. Solving the latter simultaneously with 
Eqs. HA.8I) and l|A.9li . which do not depend on P as well, 
allows to determine p(T), 0(T), and z(T) defining the 
spinodal line. Then one can determine £(T) by means of 
Eq. ES| , and finally P(T) = Tlnf(T). 
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